Human connectome project data.
Preprocessing steps involved ( at this instructional level ) for each subject:
rdir = "/Users/stnava/code/structuralFunctionalJointRegistration/"
if ( ! exists("id") ) id = '3026'
if ( ! exists("id") ) id = '2001'
# collect image data
# data/LS2001/unprocessed/3T/T1w_MPR1/LS2001_3T_T1w_MPR1_gdc.nii.gz
t1fn = paste0( rdir, 'data/LS',id, "/unprocessed/3T/T1w_MPR1/LS",id,"_3T_T1w_MPR1_gdc.nii.gz")
# now the bold data
boldfnsL = Sys.glob( paste0( rdir, 'data/LS',id, "/LS", id, "fmri/unprocessed/3T/rfMRI_REST1_*/*REST1_LR_gdc.nii.gz" ) )
boldfnsR = Sys.glob( paste0( rdir, 'data/LS',id, "/LS", id, "fmri/unprocessed/3T/rfMRI_REST1_*/*REST1_RL_gdc.nii.gz" ) )
# get the ref data
reffns = Sys.glob( paste0( rdir, 'data/LS',id, "/LS", id, "fmri/unprocessed/3T/rfMRI_REST1_*/*SBRef_gdc.nii.gz" ) )
Let us first “undistort”
if ( ! exists( "und" ) ) {
i1 = antsImageRead( reffns[1] )
i2 = antsImageRead( reffns[2] )
und = buildTemplate( i1, list( i1, i2 ), "SyN" )
}
t1 = antsImageRead( t1fn ) %>% n3BiasFieldCorrection( 8 ) %>%
n3BiasFieldCorrection( 4 )
bmask = getMask( und, lowThresh = mean( und ) * 0.75, Inf, 3 ) %>% iMath("FillHoles")
# this is a fragile approach - not really recommended - but it is quick
t1mask = getMask( t1, lowThresh = mean( t1 ) * 1.1, Inf, 5 ) %>% iMath("FillHoles")
t1rig = antsRegistration( und * bmask, t1 * t1mask, "BOLDRigid" )
t1reg = antsRegistration( und * bmask, t1 * t1mask, "SyNOnly", initialTransform = t1rig$fwd,
synMetric = 'CC', synSampling = 2, regIterations = c(5) )
###########################
Use the BOLD mask to extract the brain from the t1 (for expedience, usually we would have run antsCorticalThickness already.)
t1maskFromBold = antsApplyTransforms( t1, bmask, t1reg$invtransforms,
interpolator = 'nearestNeighbor' )
t1 = n4BiasFieldCorrection( t1, t1mask, 8 ) %>%
n4BiasFieldCorrection( t1mask, 4 )
bmask = antsApplyTransforms( und, t1mask, t1reg$fwdtransforms,
interpolator = 'nearestNeighbor' ) %>% morphology("close",3)
ofn = paste0( rdir, "features/LS", id, '_mask.nii.gz' )
antsImageWrite( bmask, ofn )
t1toBOLD = antsApplyTransforms( und, t1, t1reg$fwdtransforms )
ofn = paste0( rdir, "features/LS", id, '_t1ToBold.nii.gz' )
antsImageWrite( t1toBOLD, ofn )
plot( t1, t1mask, alpha = 0.5, axis = 3 )
## NULL
############################
here we might repeat the registration process – if we trust the masks
Exercise: Try it yourself.
Now we can segment the T1.
# a simple method
################################################
if ( ! exists( "t1seg" ) ) {
qt1 = iMath( t1, "TruncateIntensity", 0, 0.95 )
t1seg = kmeansSegmentation( qt1, 3, t1mask, 0.2 )
volumes = labelStats( t1seg$segmentation, t1seg$segmentation )
rownames( volumes ) = c("background",'csf', 'gm', 'wm' )
volumes$NormVolume = volumes$Volume / sum( volumes$Volume[-1])
pander::pander( volumes[ , c("LabelValue","Volume","NormVolume")] )
# if we look and realize this is not good - fix & redo - caused by local hyperintensity
if ( volumes["wm","NormVolume"] < 0.2 ) {
t1mask2 = t1mask * thresholdImage( t1seg$segmentation, 1, 2 )
t1seg = kmeansSegmentation( t1, 3, t1mask2, 0.1 )
}
}
plot( t1, t1seg$segmentation, axis = 3, window.overlay=c(0,3) )
## NULL
boldseg = antsApplyTransforms( und, t1seg$segmentation, t1reg$fwdtransforms,
interpolator = 'nearestNeighbor' )
plot(und,boldseg,axis=3,alpha=0.5)
## NULL
########################################
if ( ! exists( "treg" ) ) {
data( "powers_areal_mni_itk" )
myvoxes = 1:nrow( powers_areal_mni_itk )
anat = powers_areal_mni_itk$Anatomy[myvoxes]
syst = powers_areal_mni_itk$SystemName[myvoxes]
Brod = powers_areal_mni_itk$Brodmann[myvoxes]
xAAL = powers_areal_mni_itk$AAL[myvoxes]
if ( ! exists( "ch2" ) )
ch2 = antsImageRead( getANTsRData( "ch2" ) )
treg = antsRegistration( t1 * t1mask, ch2, 'SyN' )
}
concatx2 = c( treg$invtransforms, t1reg$invtransforms )
pts2bold = antsApplyTransformsToPoints( 3, powers_areal_mni_itk, concatx2,
whichtoinvert = c( T, F, T, F ) )
## Warning in data.matrix(points): NAs introduced by coercion
## Warning in data.matrix(points): NAs introduced by coercion
## Warning in data.matrix(points): NAs introduced by coercion
## Warning in data.matrix(points): NAs introduced by coercion
## Warning in data.matrix(points): NAs introduced by coercion
ptImg = makePointsImage( pts2bold, bmask, radius = 3 )
plot( und, ptImg, axis=3, window.overlay = range( ptImg ) )
## NULL
bold2ch2 = antsApplyTransforms( ch2, und, concatx2,
whichtoinvert = c( T, F, T, F ) )
# check the composite output
plot( bold2ch2 , axis=3 )
## NULL
###########################################################
back to the fmri …
* undistort
* motion correction
then we will be ready for first-level analysis …
if ( ! exists( "motcorr" ) ) {
bold = antsImageRead( boldfnsR )
avgBold = getAverageOfTimeSeries( bold )
# map to und
boldUndTX = antsRegistration( und, avgBold, "SyN" )
boldUndTS = antsApplyTransforms( und, bold, boldUndTX$fwd, imagetype = 3 )
avgBold = getAverageOfTimeSeries( boldUndTS )
motcorr = antsrMotionCalculation( boldUndTS, verbose = F, typeofTransform = 'Rigid' )
}
plot( ts( motcorr$fd$MeanDisplacement ) )
print( names( motcorr ) )
## [1] "moco_img" "moco_params" "moco_avg_img" "moco_mask"
## [5] "fd" "dvars"
trim the first \(k\) time points
trim the high motion time points and their \(k\) neighbors
goodtimes = rep( TRUE, dim( motcorr$moco_img )[4] )
goodtimes[ 1:10 ] = FALSE # signal stabilization
highMotionTimes = which( motcorr$fd$MeanDisplacement >= 0.5 )
for ( highs in -2:2 )
highMotionTimes = c( highMotionTimes, highMotionTimes + highs )
highMotionTimes = sort( unique( highMotionTimes ))
goodtimes[ highMotionTimes ] = FALSE
print( table( goodtimes ) )
## goodtimes
## FALSE TRUE
## 10 410
We can then perform regression only in the “good” time points.
Or we can impute / interpolate ( see the RestingBold article/vignette in ANTsR documentation ).
now we can do some additional level-one analysis to extract relevant networks.
* default mode
* motor
# use tissue segmentation to guide compcor
# FIXME - provide reference for this approach
csfAndWM = thresholdImage( boldseg, 1, 1 ) + ( thresholdImage( boldseg, 3, 3 ) %>% iMath("ME",1))
ccmat = timeseries2matrix( motcorr$moco_img, csfAndWM )
noiseu = compcor( ccmat, ncompcor = 10 )
# convert GM to a matrix
smth = c( 5, 5, 5, 2.0 ) # this is for sigmaInPhysicalCoordinates = T
smth = rep( 1, 4 ) # might need to mess with this some ...
simg = smoothImage( motcorr$moco_img, smth, sigmaInPhysicalCoordinates = F )
gmmat = timeseries2matrix( simg, thresholdImage( boldseg, 2, 2 ) )
tr = antsGetSpacing( bold )[4]
gmmatf = frequencyFilterfMRI( gmmat, tr = tr, freqLo = 0.01, freqHi = 0.1, opt = "trig" )
##
## 'mFilter' version: 0.1-3
##
## 'mFilter' is a package for time series filtering
##
## See 'library(help="mFilter")' for details
##
## Author: Mehmet Balcilar, mbalcilar@yahoo.com
# smoothing or filtering?
dfnpts = which( powers_areal_mni_itk$SystemName == 'Default Mode' ) # &
# powers_areal_mni_itk$Anatomy == "Posterior Cingulate Gyrus" )
dfnmask = maskImage( ptImg, ptImg, level = dfnpts, binarize = T )
dfnmat = timeseries2matrix( simg, dfnmask )
dfnmatf = frequencyFilterfMRI( dfnmat, tr = tr, freqLo = 0.01, freqHi = 0.1, opt = "trig" )
dfnsignal = rowMeans( data.matrix( dfnmatf ) )
dfndf = data.frame( signal = dfnsignal, noiseu = noiseu, fd = motcorr$fd )
myform = paste( " mat ~ ", paste( names( dfndf ), collapse = "+") )
dfnmdl = ilr( dfndf[goodtimes,],
list( mat = data.matrix( data.frame( gmmatf ))[goodtimes,] ), myform )
####################################
##############
now we can look at the full continuous beta map for the default mode network
dfnBetaImg = makeImage( thresholdImage( boldseg, 2, 2 ), dfnmdl$tValue["signal",] )
plot( und, dfnBetaImg, alpha = 1.0, axis = 3, window.overlay = c(3, max(dfnBetaImg )),
nslices = 24, ncolumns = 6 )
## NULL
ofn = paste0( rdir, "features/LS", id, '_dfnBetaImg.nii.gz' )
antsImageWrite( dfnBetaImg, ofn )
ofn = paste0( rdir, "features/LS", id, '_undistort.nii.gz' )
antsImageWrite( und, ofn )
ofn = paste0( rdir, "features/LS", id, '_t1ToBold.nii.gz' )
antsImageWrite( t1toBOLD, ofn )
ofn = paste0( rdir, "features/LS", id, '_mask.nii.gz' )
antsImageWrite( bmask, ofn )
Now repeat all of this for the next subject by changing the ID.
Finally, use all of this data to do a joint registration using both structural and functional features.
# antsRegistration with multivariateExtras
id1 = '2001'
s1f1 = antsImageRead( paste0( rdir, "features/LS", id1, '_dfnBetaImg.nii.gz' ) )
s1f2 = antsImageRead( paste0( rdir, "features/LS", id1, '_undistort.nii.gz' ) )
s1f3 = antsImageRead( paste0( rdir, "features/LS", id1, '_t1ToBold.nii.gz' ) )
s1mask = antsImageRead( paste0( rdir, "features/LS", id1, '_mask.nii.gz' ) )
id2 = '3026'
s2f1 = antsImageRead( paste0( rdir, "features/LS", id2, '_dfnBetaImg.nii.gz' ) )
s2f2 = antsImageRead( paste0( rdir, "features/LS", id2, '_undistort.nii.gz' ) )
s2f3 = antsImageRead( paste0( rdir, "features/LS", id2, '_t1ToBold.nii.gz' ) )
s2mask = antsImageRead( paste0( rdir, "features/LS", id2, '_mask.nii.gz' ) )
#############
jrig = antsRegistration( s1f3 * s1mask, s2f3 * s2mask, "Affine" )
jreg = antsRegistration( s1f3, s2f3, "SyNOnly", initialTransform = jrig$fwd, mask = s1mask,
multivariateExtras = list(
list( "mattes", s1f2, s2f2, 1, 32 ),
list( "mattes", s1f1, s2f1, 1, 32 ) ), verbose = FALSE )
#############